function [ nu, taw ] = ep_params ( K, y )
%This funciton will estimate the parameters of the ep for Gaussian Process
%efined by the kernel matrix K and the predictions vector y

n = size(K, 1);
taw = zeros(n, 1);
nu =  zeros(n, 1);
mu =  zeros(n, 1);
Sig = K;

converged = false;
while ~ converged
    prev_nu = nu;
    prev_taw = taw;
    for i = 1 : n
        taw_i = Sig(i,i) - taw(i);
        nu_i  = Sig(i,i) * mu(i) - nu(i);
        if ((Sig(i,i) - taw_i) == 0)
           sig_i = 1000;
        else            
            sig_i = 1. / (Sig(i,i) - taw_i);
        end
        mu_i = sig_i * (Sig(i,i) * mu(i) - taw(i) * nu(i));
        zi = y(i) * mu_i / sqrt(1 + sig_i);
        cdfz = normcdf(zi);
        if (cdfz == 0)
            cdfz = 1e-4;
        end
        mu_hat_i    = mu_i + y(i) * sig_i * normpdf(zi) / (cdfz * sqrt(1 + sig_i) );
        sigma_hat_i = sig_i + (sig_i * sig_i * normpdf(zi) / (1 + sig_i) * cdfz) * (zi + normpdf(zi)/cdfz);
        delta_taw = sigma_hat_i - taw_i - taw(i);
        taw(i) = taw(i) + delta_taw;
        nu(i) = sigma_hat_i * mu_hat_i - nu_i;
        Sig = Sig - 1./(1./delta_taw + Sig(i,i)) * Sig(:,i) * Sig(:,i)';
        mu = Sig * nu;
    end
    S = diag(taw);
    S = sqrtm(S);
    L = chol(eye(n,n) + S * K * S,'lower');
    V = L'\S*K;
    Sig = K - V'*V;
    mu = Sig * nu;
    
    if ( norm(nu - prev_nu) < .1 & norm(taw - prev_taw) < .1 )
        converged = true;        
    end 
end 


end

